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ABSTRACT 


This  report  deals  with  minimax  approximations  to  functions  defined  on  the  real  line.  These 
approximations  are  of  particular  use  to  mathematicians  having  access  to  a  digital  computer 
facility. 

The  cases  considered  are  for  functions  defined  as  discrete  points  on  the  Cartesian  plane  and 
for  continuous,  bounded  functions  defined  on  a  closed  interval  of  the  real  line.  A  method  for 
obtaining  the  approximations  using  linear  programming  techniques  is  presented  for  each  case. 

An  analysis  of  the  error  function  and  a  brief  outline  of  the  computer  program  are  included. 
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SECTION  I 

INTRODUCTION 


The  general  problem  considered  in  this  paper  is  the  approximation  of  a  function  f  by  an 
approximating  function  F.  The  re  are  two  main  reasons  for  being  interested  in  approximating 
a  given  function.  In  the  first  place,  if  we  are  using  a  digital  computer,  nonelementary  functions 
may  take  considerable  machine  time  to  compute,  and  if  we  are  to  evaluate  such  a  function  a 
large  number  of  times,  it  will  be  advantageous  to  replace  it  with  a  more  elementary  one.  In  the 
second  place,  if  the  function  f  is  given  in  tabular  form,  a  good  approximating  function  will  re¬ 
place  the  need  for  table  look-up  and  interpolation  and  will  economize  on  storage. 

In  approximating  a  given  function  we  are  faced  with  two  major  problems:  (1)  to  determii  e 
the  form  of  the  approximating  function  and  (2)  to  determine  a  criterion  for  the  accuracy  of  the 
approximation.  The  primary  consideration  in  determining  the  type  of  approximating  function 
to  use  is  that  this  function  be  compatible  with  the  function  to  be  approximated.  Thus  we  would 
not  approximate  a  function  which  has  a  singularity  with  a  polynomial.  There  may  be,  however, 
certain  practical  considerations  which  influence  the  choice  of  the  approximating  function.  In 
the  following  we  limit  ourselves  to  polynomials  and  functions  which  we  can  approximate  with 
polynomials,  since  they  are  easily  evaluated. 

Just  as  there  are  numerous  choices  for  the  approximating  function,  there  are  many  ways  of 
measuring  the  accuracy  of  the  approximation.  Among  the  more  common  methods  are  exact 
matching  at  the  sample  points;  least- squares,  which  is  well  suited  to  data  subject  to  a  random 
error;  and  Chebycheff,  which  reduces  the  maximum  error  to  a  minimum.  The  actual  choice  of 
the  criterion  to  use  depends  on  the  particular  problem.  The  Chebycheff  criterion  is  partic¬ 
ularly  useful  when  the  noise  connected  with  the  data  is  negligible,  or  when  it  is  desired  to 
keep  the  maximum  error  minimal.  We  will  concern  ourselves  with  the  Chebycheff  criterion. 

The  problem  we  will  consider  is  that  of  finding,  among  all  polynomials  Q  of  degree  less  than 
or  equal  to  n,  a  polynomial  P  such  that  given  a  function  f  defined  on  a  set  S  of  real  numbers 

max  |  f(x)  -  P(x)  |  =  min  (max  |f(x)  -  Q(x)  |  ). 

xcS  Q  xcS 

The  set  S  may  be  a  discrete  set  of  points  or  an  interval  on  the  real  line.  P  is  said  to  b-^  the 
polynomial  which  best  approximates  f  on  S  in  the  Chebycheff  sense. 

We  will  show  the  existence  and  uniqueness  of  such  a  polynomial  P,  describe  the  character¬ 
istics  of  the  error  function,  and  develop  the  necessary  computational  methods  to  find  P. 
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SECTION  II 

CHEBYCHEFF  APPROXIMATION  TO  A  FUNCTION 
DEFINED  ON  A  FINITE  POINT  SET 


Given  a  function  f  defined  on  a  point  set  S  of  m  points  of  the  real  line,  wc  want  to  find  a  poly¬ 
nomial  P,  of  degree  less  than  or  equal  to  n,  n+2<m,  for  which  the  error 

max  |  f(x)  -  P(x)  | 
xeS 

is  a  minimum.  We  call  P  the  polynomial  of  degree  less  than  or  equal  to  n  of  best  approximation 
to  f  on  S. 

The  following  theorem  of  C.  de  la  Vallee  Poussin  (Reference  19)  characterizes  the  poly¬ 
nomial  P. 

Theorem:  Let  P  be  a  polynomial  of  degree  less  than  or  equal  to  n,  of  best  approximation  to 
a  function  f  defined  on  a  point  set  S  of  m>n+2  points. 

Then  P  is  uniquely  characterized  by  the  existence  of  at  least  n+2  points  of  S  at  which  the 
difference  |f(x)  -  P(x)|  assumes  its  maximum  value  on  S,  with  the  algebraic  sign  of  f(x)  -  P(x) 
alternating  at  these  points. 

Proof:  We  suppose  |  f(x)  -  P(x)  |  assumes  its  maximum  value  at  no  more  than  n+1  points  of 
the  set  S  and  arrive  at  a  contradiction. 

Let  S/  be  this  set  of  at  most  n+1  points.  We  can  then  find  a  polynomial  Q  of  degree  less  than 
or  equal  to  n  such  that  Q(x)  has  the  same  algebraic  sign  as  f(x)  -  P(x)  for  every  point  x  in  S  . 

Let  E(P;S)  =  max  I  f(x)  -  P(x)  I 
xcS 

E(P;S-S')  =  max|f(x)  -  P(x)| 
xcS-S' 

and  II  Q(S)  ||  =  max|Q(x)|  . 

X€S 

Since  E(P;S)>E(P;S-S')  we  can  find  anc>0  such  that  cl lQ(S) 1 1  <  E(P;S)  -  E(P;S-S'). 

Let  R(x)  =  P(x)  +  eQ(x).  R(x)  is  then  a  polynomial  of  degree  less  than  or  equal  to  n.  We  now 
show  that  under  the  hypothesis  that  there  are  at  most  nfl  points  at  which|f(x)  -  P(x)|assumes 
its  maximum  value  on  S,  the  polynomial  R(x)is  a  better  approximation  to  f(x)  than  P(x),  which 
is  a  contradiction  since  P(x)  was  taken  as  the  best  approximation  to  f(x)  on  S. 

On  the  subset  S1  of  S  we  have 

I  f(x)  -  R(x)  I  =  I  f(x)  -  P(x)  -«Q(x)  I  =  I  f(x)  -  P(x)  I  -  €  I  Q(x)  I. 

The  last  equality  holds  because  Q(x)  has  the  same  algebraic  sign  as  f(x)  -  P(x)  on  s'. 

But  |  f (x)  -  P(x)  |  -  e  |Q(x)  |  <  |  f(x)  -  P(x)|  .  Hence  for  xeS*  we  have  |f(x)  -  R(x)  |  < 

|f(x)  -  P(x)|  . 
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On  the  subset  S-S*  of  S  we  have 

|  f(x)  -  R(x)  |  =  |  f(x)  -  P(x)  -«Q(x)l<  lf(x)  -  P(x)  I  +€llQ(S)  II  . 
Then 

max  I  f(x)  -  R(x)  I  <  max  I  f(x)  -  P(x)  I  +  dlQ(S)  II  . 
x«S-S'  xcS-S1 

But 

max  I  f(x)  -  P(x)  I  +  c  II  Q(S)  II  <  E(P;S-s')  +  (E(P;S)  -  E(P;S-S‘)  ) 
xeS-S1 


<E(P;S) 

Hence  max  I  f(x)  -  R(x)  l<E(P;S). 

X€S-S' 

Since  S  is  the  union  of  S'  and  (S-S1)  the  polynomial  R(x)  gives  a  better  approximation  to  f(x) 
on  S  than  P(x),  which  contradicts  the  hypothesis  that  P(x)  is  the  polynomial  of  best  approxi¬ 
mation. 

Having  characterized  the  polynomial  P  of  best  approximation,  we  now  show  it  is  unique. 

Theorem:  Let  P  be  a  polynomial  of  degree  less  than  or  equal  to  n,  of  best  approximation 
to  a  function  f  defined  on  a  set  S  of  m>n+2  points.  Then  P  is  unique. 

Proof:  Assume  the  contrary,  that  there  exists  a  polynomial  Q  of  degree  less  than  or  equal 
to  n,  different  from  P,  such  that  max  I  f(x)  -  P(x)  I  =  max  I f(x)  -  Q(x)  I. 

xcS  xcS 

Let  E(S)  =  max  I  f(x)  -  F(x)  I  =  max  I  f(x)  -  Q(x)  I . 
xcS  xeS 


Then  for  xcS 

-E(S)<  f (x)  -  P(x)  <  E(S) 
and  -E(S)  £  f(x)  -  Q(x)  1  E(S). 

Addition  and  division  by  two  yields 

-E(S)  <f(x)  -  .5(P(x)  +  Q(x)  )  <E(S). 

Then  R(x)  =  .5(P(x)  -  Q(x)  )  is  also  a  polynomial  of  degree  less  than  or  equal  to  n  of  best 
approximation  to  f(x)  on  S. 

Consider  a  point  x1  of  S  for  which  f(x')  -  L(x')  =  E(S), 

Then  2(f(x')  -  R(x')  )  =  (f(x')  -  P(x')  )  +  (f (x1 )  -  Q(x')  )  =  2E(S). 

Since  both  f(x')  -  P(x')  and  f(x' )  -  Q(x')  are  less  than  or  equal  to  E(S)  we  must  have 

f(x‘)  -  P(x')  =  f(x')  -  Q(x'),  i.e.,  P(x')  =  Q(x'). 
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It  can  likewise  be  shown  that  for  any  point  x"  of  S  for  which  f(x")  -  R(x")  =  -E(S)  we  will 
have  P(x")  =  Q(x"). 

By  the  previous  theorem  there  are  at  least  n+2  points  at  which  f(x)  -  R(x)  =  ±E(S).  Hence 
P(x)  and  Q(x)  coincide  at  a  minimum  of  m+2  points  of  S.  But  then  P(x)  and  Q(x)  are  identical 
since  two  polynomials  of  degree  less  than  or  equal  to  n,  coinciding  at  n+2  points,  coincide 
everywhere. 

We  now  return  to  the  problem  of  finding  the  unique  polynomial  P,  of  degree  less  than  or 
equal  to  n,  of  best  approximation  to  a  function  f  defined  on  a  set  S  of  m  points,  where  m  is 
greater  than  or  equal  to  n+2. 

The  problem  is  one  of  determining  the  coefficients  c„,  .  .  .  ,  c  and  the  deviation  E(S)  such 
that  U  n 


j=n 

t  (S)  =  max  I  f  (x  )  -  X  c  •  x  M 
X€  S  j  =  0  1 


is  a  minimum. 

=  J 


Let  a.j  =  x^,  l<ilm  and  01  jin.  The  problem  is  then  to  find  cQ,  ^ . cr  and  E(S)  such  that 


E(S)  +  Z  a.,  c.  >  f 
j=0  'J  1  ' 


i  =  I, 


m 


E(S)  -  X  a*,  c.  >  -f. 
j  =  0  'I  i  1 


E(S)  >  0 


and  E(S)  is  a  minimum. 

We  now  have  a  system  of  2m+l  linear  inequalities  in  the  n+2  unknowns  c^,  c^ . and 

E(S) ,  and  a  linear  functional  Z  =  E(S)  to  be  minimized. 

We  will  see  that  there  is  an  efficient  way  of  adapting  our  approximation  problem  to  a  linear 
programming  form.  To  show  this  we  define  certain  terms  and  concepts  of  linear  programming. 

Associated  with  every  linear  programming  problem  is  another  linear  programming  problem 
called  the  ‘dual’.  Below  is  given  the  relation  between  a  typical  ‘primal’  problem  and  its  dual. 
It  should  be  noted  that  either  of  the  two  problems  may  be  called  the  primal,  and  the  other  would 
then  be  the  dual. 


Primal  Problem:  Find  a  column  vector  X  =  (x.,  . 
functional 


x  )  which  maximizes  the  linear 

n  ' 


z  =  cx 


subject  to  the  restrictions 

AX  =  B  and  X>0 
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where  C  =  (c^,  ....  c)  is  a  row  vector,  B  =  (bj, 
is  a  p  by  q  matrix. 


b  )  is  a  column  vector,  and  A  =  (a..) 
p  ij 


Dual  Problem;  Find  a  row  vector  Y  =  (y^, 

W  =  YB 


y  )  which  minimizes  the  linear  functional 


subject  to  the  restriction  YA  ^  C,  where  A,  B,  and  C  are  defined  as  in  the  primal  problem. 
Note  that  there  is  no  restriction  on  the  sign  of  Y. 

The  following  is  the  well-known  duality  theorem,  which  establishes  the  important  relation 
between  a  primal  problem  and  its  dual. 

We  define  a  feasible  solution  as  one  which  satisfies  the  constraints  of  the  given  problem,  and 
an  optimal  solution  as  a  feasible  solution  which  maximizes  or  minimizes  the  object  function  Z 
or  W,  as  the  case  may  be. 

Theorem:  If  either  the  dual  or  the  primal  problem  has  a  finite  optimal  solution,  the  other 
problem  has  a  finite  optimal  solution,  and  the  extremes  of  the  linear  functionals  are  equal, 
i.e.,  max  Z  =  min  W. 


A  proof  of  this  theorem  may  be  found  in  thj  book,  Linear  Programming,  by  S.  I.  Gass 

(Reference  7). 


We  solve  the  above  linear  programming  problem  using  the  Simplex  Method  due  to  G.  B. 
Xnt  ig  (Reference  5).  A  ‘brsic’  feasible  solution  is  one  with,  at  most,  p  positive  variables, 

,v:!  o;ie  which  expresses  B  in  terms  ofp  linearly  independent  vectors  P^,  j  =  1 . p,  where 

a  .)  is  the  .th  column  of  A. 


(a 


j’ 


ir 


p  r 


J 


The  algorithm  starts  with  a  basic  feasible  solution.  We  can  always  find  such  a  solution, 
M.ough  it  may  be  necessary  to  introduce  certain  artificial  vaiiables  to  do  so.  If  this  solution 
not  optimal,  we  can  either  find  a  class  of  solutions  in  which  there  is  no  bound  on  the  object 
emotion  or  we  can  construct  a  new  basic  feasible  solution  in  which  the  value  of  the  object 
.unction  is  not  less  than  that  for  the  previous  solution.  The  simplex  algorithm  consists  of  re¬ 
peating  this  cycle  until  we  have  constructed  an  optimal  basic  feasible  solution  or  a  class  of 
feasible  solutions  for  which  the  object  function  is  unbounded. 


Let  G  be  the  p  by  p  matrix  formed  by  the  p  linearly  independent  vectors  P.  from  A  which 

i  ^ 

correspond  to  the  optimal  feasible  solution.  We  call  G  the  optimal  basis.  Let  X  be  the  p  dimen¬ 
sional  optimal  basic  feasible  solution.  At  the  termination  of  the  simplex  algorithm  there  will 
be  a  p  dimensional  row  vector  C'  such  that 

X'  =  G_1B  and  max  Z  =  C'X’ 


If  we  let  Y'  be  the  optimal  basic  feasible  solution  of  the  dual,  then 
min  W  =  Y'B  =  C'G_1B  =  C'X'  =  max  Z,  i.e.,  Y'  =  C'G-1. 

Most  standard  linear  programming  codes  for  digital  computers  will  yield  the  solution  to  both 
the  primal  and  the  dual  problems. 
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Consider  now 


our 


approximation  problem. 


Since  a. 

10 


o 

x. 

1 


we  let  the  matrix  A  be 


1 1  > 


i  n  » 


1 

1 

1 

> 

I 

1 

) 

» 

1 

- 1 

°2. 

1 

ami  * 

_Q.  . 

a2n 

1 

» 

a  mn  > 

c 

.  a 
.  i 

)  ■ 

,0 

,0 


-a 


mn 


,0 


We  also  let 


and 


F  =  (f. . f  ,  -  f  ,  . 

1  m  1 


Y  =  (E(S) ,  Cq . cn) 


0) 


wh^re  U  is  a  2m+l  by  1  column  vector.  The  problem  is  then  to  find  a  row  vector  Y  which  mini¬ 
mizes  the  linear  functional  W  =  YU  subject  to  YA>F,  which  is  precisely  the  dual  form  de¬ 
scribed  above.  We  can  consider  the  primal  problem  as  one  of  finding  a  column  vector  X  which 
maximizes  Z  =  FX  subject  to  AX  =  U  and  X>0. 

We  can  now  solve  this  primal  version  of  the  problem  using  a  standard  linear  programming 
code,  which  will  yield  a  solution  to  the  dual  as  well. 

There  are  three  reasons  for  taking  the  “dual  approach”  rather  than  solving  for  Y  directly: 

1.  Since  Y  is  not  restricted  in  sign,  we  would  be  forced  to  add  variables  to  our  problem 
before  attempting  a  solution.  The  dual  approach  eliminates  this  need. 

2.  The  simplex  method  is  most  efficient  in  terms  of  time  when  the  number  of  variables  is 
greater  than  the  number  of  equations.  Since  in  our  problem  we  have  n+2  variables  but  2m+l 
equations,  where  m>n+2,  it  is  advantageous  to  consider  the  dual  so  that  we  have  2m+l  variables 
and  n+2  equations. 

3.  The  addition  of  points  at  which  we  want  to  approximate  f(x)  corresponds  to  the  addition 
of  variables  in  the  problem  as  we  consider  it,  and  this  can  be  done  without  restarting  the 
algorithm  at  its  initial  phase.  But  addition  of  points  would  mean  the  addition  of  constraints  if 
we  solve  our  problem  directly.  Having  obtained  a  solution  vector  which  satisfies  the  original 
constraints,  we  would  have  to  start  the  solution  anew  to  take  into  consideration  the  new  con¬ 
straints. 


Interpolation  Error:  Let  f  be  a  continuous  function  with  at  least  n+2  derivatives,  defined 
on  a- iet~S~oT~m  points.  Having  obtained  a  polynomial  P,  of  degree  less  than  or  equal  to  n, 
m>n+2,  which  is  a  best  approximation  to  f(x),  we  may  want  to  use  P  as  an  interpolating  poly¬ 
nomial  to  interpolate  between  the  given  points  of  S. 
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We  will  call  the  points  of  S  at  which  |  f(x)  -  P(x)  |  takes  on  its  maximum  value  the  critical 
points,  and  we  assume  that  there  are  exactly  n+ 2  of  them.  We  further  assume  that  the  smallest 
and  largest  of  the  points  in  S  are  in  the  critical  set.  Let  x'  and  x"  be  these  two  points.  Given 
these  conditions,  we  can  obtain  an  interpolation  error  for  P  over  [x1,  x"] . 

Let  e(i)  =  f(x.)  -  P(x.),  i  =  1 . n+2,  for  x.  in  the  critical  set.  We  can  construct  a  poly¬ 

nomial  L(x)  of  degree  n+1  which  will  assume  the  values  e(i)  at  the  points  x.  of  the  critical  set. 
We  use  the  Lagrangian  interpolation  formula  to  obtain 

L(x)  =  h(l,  x)  e(l)  +  .  .  .  +  h(n+2,  x)  e(n+2) 
where 


h  ( i ,  x  )  = 


(X  -  Xj)  .  .  .  (X  -  X.  _  j)  (X  -  X.  +  j)  .  .  .  (X  -  xn+2) 

(X.  -  X.)  .  .  .  (X.  -  X.  .)  (X.  -  X.  -  X.  .}  ...  (X.  -  X  _) 

v  i  r  i  i  -  r  v  i  i  i+l'  v  i  n  +  2 


The  error  term  of  this  formula  is 


RU  ) 


G(n+2)(v)  Q(x) 
(n+2) 


where  Q(x)  =  (x  -  x^  .  .  .  (x  -  *n+2)- 


The  x^’s  are  the  critical  points,  and  G(x)  =  f(x)  -  P(x),  and  the  poin*  v  is  some  point  in 

[x’f  x"] .  A  complete  discussion  of  the  Lagrangian  interpolation  formula  and  the  error  term 
can  be  found  in  (Reference  8). 

Since  P  is  a  polynomial  of  degree  at  most  n,  its  (n+2)nd  derivative  is  identically  zero;  hence, 


G(n+2)(V)  =  f(n+^(V)  and  R(x)  = 


f(n+2)(v)  Q(x) 
(n+2)! 


The  error  f(x)  -  P(x)  is  then  equal  to  L(x)  +  R(x),  and  this  is  the  error  incurred  when  using 
P(x)  as  an  interpolating  polynomial. 
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SECTION  m 

CHEBYCHEFF  APPROXIMATION  TO  A  FUNCTION 
DEFINED  ON  A  CLOSED  INTERVAL 
ON  THE  REAL  LINE 


Let  f  be  a  given  function  defined  and  continuous  on  a  closed  interval  [  a,  b].  We  are  concerned 
with  finding  a  polynomial  P,  of  degree  less  than  or  equal  to  n,  such  that  the  error 

max  |  f(x)  -  P(x)  I 
a<x<b 


is  a  minimum. 

The  following  theorem  of  Chebycheff  (Reference  1)  characterizes  the  polynomial  P. 

Theorem:  Among  all  the  polynomials  of  degree  less  than  or  equal  to  n  defined  on  the  real 
line,  the  polynomial  P  which  deviates  least  from  a  given  continuous  function  is  uniquely  char¬ 
acterized  by  the  fact  that  in  the  interval  [a,  b  ]  over  which  we  approximate  the  given  function, 
the  number  of  consecutive  points  at  which  the  difference  f(x)  -  P(x),  with  alternating  signs, 
assumes  the  value 


max  |  f(x)  -  P(x)  | 
a<x<b 

is  not  less  than  n+2. 

The  next  theorem  is  important  for  computational  reasons.  We  define  a  local  maximum  as 
follows:  we  say  f  has  a  local  maximum  at  a  point  x  if  there  exists  a  8>0  such  that  f(y)  <  f(x) 
for  all  y  such  that  0  <  I  x  -  y  I  <  B  .  A  local  minimum  can  be  defined  similarly  by  reversing 
the  inequality  between  f(x)  and  f(y).  We  now  show  that  the  function 

E(P;  [a,  b]  )  =  I  max  f(x)  -  P(x)  I  , 
a<x<b 

where  P(x)  =  c^  +  c^  x  +  .  .  .  =  cnxn  is  an  arbitrary  polynomial  of  degree  ^  n,  is  convex,  and 

hence  a  local  maximum  or  minimum  is  the  unique  maximum  or  minimum  for  the  whole  interval 
[a,  b]. 

Theorem:  Let  9  be  in  [o,  l] ,  and  let  P,  and  P  ,  be  polynomials  of  degree  less  than  or  equal 
to  n.  Then 

E(0Pj  +  (1  -9)  P2;  [a,  b]  )*0E(Pi;  [a,  b]  )  +  (1  -9)  E(P2;  [a,  b]  ) 
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Proof: 


E(^P,+(  1—^)  P_  ;  [  a ,  b  1 )  =  max  I  f  (x )  -  ]T  [6c  +  (  I  -  #)  c_  ) x'  I 

1  J  a<  x<  b  i=0  i  i 

n  n  . 

=  maxi  0(  f  (  x) c,  x  1  )  +(  I -6)  ( f  (  x)  c2  x  )l 
a<x<b  i  =  0  '  i  =0  * 

n  i  n  i 

<  max  (#lf(x)-£  c,  x  I  +  (  I  ~  6)  I  f  (x )  -  £  c2  x  I) 

a<x<  b  i=0  i  i=0  i 

n  n 

<9  maxi  f  (  x)  c(  x'  1+  (  1-0)  max  if  ( x )  —  c2  x'  I 
a<x<b  •  =0  '  a<  x<  b  *=0  1 


=  6 E(  P,  ;  [o,b] )  +( l-0)E(P2  ;  [a  ,b]  ). 


Let  us  suppose  now  that  the  function  E  has  local  minima  for  polynomials  Pj  and  of  degree 
less  than  or  equal  to  n,  and  let  Q  be  in  [o,  1  ].  We  have  just  proved  that 

E(0Pj  +  (1  -6)  P2;  [a,  b]  )  <  0E(Pl5  [a,  b]  )  +  (1  -6)  E(P2;  [a,  b]  ). 

Without  loss  of  generality,  suppose  E(P^  [a,  b]  )  <  E(P2;  [a,  b]  ).  Then  by  letting  Q  approach 
0  we  can  find  a  P3  as  near  to  P2  as  we  like  such  that  E(P3;  [a,  b]  )  <  E(P2;  [a,  b]  ),  con¬ 
tradicting  the  hypothesis  that  E  has  a  local  minimum  to  P2<  By  P^  being  near  to  P2  we  mean 


(  Y.  I  c3  -  c2  )2)  l/2<  c  ,  i  =  0 
i=0  '  « 


n  , 


for  arbitrary  c>0. 

On  the  other  hand,  if  E(P^;  [a,  b]  )  =  E(P2;  [a,  b]  ),  then  there  can  be  no  P3  such  that 


C1  <  c3  <  c2  or  ci  >  c3  >  c2  »  i  =  0,  .  .  .  ,  n, 

l  1  i  i  l  i  l 

and  such  that  E(P3;  [a,  b]  )  >  E(P1;  [a,  b]  )  and  E(P3;  [  a,  b  ]  )  >E(P2;  [a,  b]  ).  We  can  see 
this  by  observing  that  E(@P1  +  (1 -0)  P2;  [a,  bj  )^E(P1;|a,  b]  ). 

We  see  then  that  the  function  E  has  no  more  than  one  minimal  value;  and  hence,  once  we 
have  a  minimum  for  E,  we  know  this  is  the  unique  minimal  value  of  the  function. 
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We  now  develop  an  algorithm  for  determining  approximately  the  unique  polynomial  P  of 
degree  less  than  or  equal  tc  n  which  is  a  best  approximation  to  a  given  function  f,  defined  and 
continuous  on  an  interval  [  a,  b].  Without  loss  of  generality,  we  can  make  the  linear  'sub¬ 
stitution  x'  =  (2x  -  a  -  b)  /  (b  -  a)  which  transforms  the  interval  [a,  b  ]  into  the  interval 

[-1,1]. 


The  exact  solution  to  the  problem  involves  the  solution  of  a  system  of  nonlinear  equations, 
and  is  quite  difficult  to  arrive  at  in  all  but  the  simplest  cases.  We  can,  however,  arrive  at  an 
approximate  solution  by  successive  approximations  on  discrete  subsets  of  the  interval  [-  1, 1  j. 

We  generate  a  sequence  Q(k;  x)  of  polynomials  of  degree  less  than  or  equal  to  n,  along  with 
a  sequence 

X(k)  =  (x(k,  i):  i  =  1,  2,  .  .  .  ,  n+2) 

of  sets  of  n+2  points  of  the  interval  [-  1,  l].  The  polynomial  Q(k;  x)  is  the  best  approximation 
to  f(x)  on  X(k),  k  =  0,  1,  2,  .  .  . 

Let 

E(Q(k;  x);  X(k))  =  max  I  Q(k;  x)  -  f(x)  |  . 
x«X(k) 

Any  set  X(0)  of  n+2  points  from  [  -  1,  l]  for  which  E(Q(0;  x);  X(0))  is  not  equal  to  zero  can 
be  chosen  to  start  the  calculations,  but  in  practice  it  has  been  found  that  faster  convergence 
is  attained  if  we  choose  the  n+2  points  at  which  the  Chebycheff  polynomial  T(n+1)  of  degree  n+1 
attains  its  extremal  values. 

Having  chosen  X(0),  we  find  Q(0;  x)  by  the  method  described  in  Section  II.  The  iterative 
method  proceeds  as  follows: 

1.  Given  Q(k;  x),  choose 

X(k  +  1)  =  (x(k  +  1,  i):  x(k  +  1,  i)  <  x(k  +  1,  i  +  1)),  i  =  0 . n+1,  so  thatQ(k;  x)  -  f(x)  has 

opposite  signs  at  adjacent  points  of  X(k+1)  and  Q(k;x)-f(x)  assumes  its  extrema  on  X(k  +  1). 
It  should  be  noted  that  X(k  +  1)  andX(k)  need  not  be  disjoint.  It  will  be  seen  that  for  most  cases 
the  two  sets  do  indeed  have  elements  in  common. 

Let 

M(k)  =  max  I  Q(k;  x)  -  f(x)  I 
x  cX{  k  +  1) 

and 

m(k)  =  min  I  Q(k;  x)  -  f(x)  I 
xcX(k  +  1) 

We  shall  prove  that  m(k)  will  be  at  least  equal  to  E(Q(k;  x);  X(k). 

2.  Q(k  +  1;  x)  and  E(Q(k  +  1;  x);  X(k  +  1)  are  then  determined  so  that  Q(k  +  1,  x)  is  the  poly¬ 

nomial  of  degree  less  than  or  equal  to  n  of  best  approximation  to  f(x)  on  X(k  +  1). 

We  accomplish  Steps  1  and  2  in  the  following  manner.  Since  the  error  Q(k;  x)  -  f(x)  takes 

on  the  values  E(Q(k;  x);  X(k))  and  -E(Q(k;  x);  X(k))  at  no  less  than  n+2  points  of  X(k),  there 
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will  be  at  least  n  +  1  points  of  [-1,  ljfor  which  Q(k;  x)  -  f(x)  is  equal  to  zero,  say 

x,'  <x  '  .  .  .  <x‘  Let  x  '  =  -1  and  x'  =  1 

1  2  n+1  0  n+2 


where  x^  <  x^ 


and  x‘  . 
n+1 


<  x 


n+2' 


The  set  X(k  +  1)  is  then  chosen  so  that  x(k  +  1,  i)  is  the  point  on  [x.’  ^,  x!J  at  which 

I  Q(k;  x)  -  f(x)  I  is  maximal.  The  set  X(k  +  1)  so  chosen  satisfies  the  requirements  in  Step  1. 
We  now  compute  Q(k  +  1;  x)  and  E(Q(k  +  1;  x);  X(k  +  1)). 


Let 


d(i)  =  I  Q(k;  x(k  +  1,  i)  -  f(x(k  +  1,  i)  I  ,  i  =  1 . n+2. 


Then 

(Q(k;  x(k  +  1,  i))  -  f(x(k  +  1;  i)))  -  (Q(k  +  1;  x(k  +  1,  i))  -  f(x(k  +  1,  i)))  = 

Q(k;  x(k  +  1,  i))  -  Q(k  +  1;  x(k  +  1,  i))  =  ±(d(i)  -E(Q(k  +  1;  x(k  +  1,  i)))). 

We  see  that  E(Q(k  +  1;  x);  X(k  +  1))  must  be  greater  than  m(k),  for  if  m(k)>E(Q(k  +  1;  x);  X(k  +  1)) 
the  polynomial  Q(k;  x)  -  Q(k  +  1;  x)  would  have  n+2  zeroes;  but  since  Q(k;  x)  and  Q(k  +  1;  x) 
are  both  of  degree  less  than  or  equal  to  n,  this  would  imply  that  Q(k;  x)  =  Q(k  +  1;  x). 


We  can  iterate  the  procedure  so  that  we  obtain 

E(Q(k;  x);  X(k))  <  m(k)  <  E(Q(k  +  1;  x);  X(k  +  1)),  k  =  0,  1,  2 . 

Let  P(x)  be  the  polynomial  of  best  approximation  to  f(x)  on  [-  1,  l] ,  and  let 

E(P;  [-  1,  l])  =  max  |  f(x)  -  P(x)  | 

-  1<X<1 


For  Q(k  +  1;  x)  not  identical  to  P(x),  the  following  relation  holds: 

E(Q(k  +1);  X(k  +  1))  <  E(P;  [-  1,  l]  )<  M(k) 

Since  X(k  +  1)  is  a  subset  of  [  -  1,  l]  ,  and  Q(k  +  1)  is  the  polynomial  of  best  approximation 
to  f(x)  on  X(k  +  1),  we  cannot  have  E(Q(k  +  1);  X(k  +  1))  >E(P;  [-  1,  1  ]  ),  for  this  would  mean 
P(x)  is  a  better  approximation  to  f(x)  on  X(k  +  1),  which  is  a  contradiction. 

It  is  also  seen  that  we  cannot  have  M(k)  <  E(P;  [  -  1,  l]  )  since  if  we  did  then  Q(k  +  1)  would 
be  a  better  approximation  to  f(x)  on  1,  1  j  than  P(x),  which  is  again  a  contradiction. 

We  now  show  that  the  process  described  above  will  indeed  result  in  convergence  of  the 
sequence  E(Q(k;  x);  X(k)),  k  =  0,  1,  2,  ...  to  the  number  E(P;  1,  lj  ). 

We  can  find  a  X  such  that  0<X<1  and 


1.  E(Q(k  +  1;  x);  X(k  +  1))  -  m(k)  >(1  -  X)  (M(k)  -  m(k)) 

Since  M(k)  >E(P;  [  -  1,  1  )  >  m(k),  we  also  have 

2.  E(Q(k  +  1;  x);  X(k  +  1))  -  m(k)  >  (1  -  X  )  (E{P;  [-  1, 

and  m(k  +  1)  >  E(Q(k  +  1;  x);  X(k  +  1))  further  implies 


)  -  m(k)) 
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3.  m(k  +  1)  -m(k)  >  (1  -  X )  (E(P;  [ -  1,  1  ] )  -  m(k)) 

Subtracting  both  sides  of  (3)  from  E(P;  [-  1,  l]  )  -  m(k),  we  obtain 

4.  E(P;  [-  1,  l]  )  -  m(k  +  1)  <  X  (E(P;  [ -  1,  1  ]  )  -  m(k)) 

But  4.  implies  that  the  quantity  E(P;  [-  1,  1  ]  )  -  m(k),  k  =  0,  1,  2 . will  tend  to  zero 

at  least  as  fast  as  the  terms  of  a  diminishing  geometric  progression  upon  repeated  applications 
of  the  algorithm.  Since  E(Q(k  +  1;  x);  X(k  1))  >  m(k),  the  quantity  E(P;  f-  1,  1  j  )  -E(Q(k  +  1;  x); 
X(k  +1))  will  likewise  tend  to  zero  at  the  same  rate. 

Note  also  that  by  1.  we  have 


M(k)  -  E(P;  [-  1,  l]  )  M(k)  -  m(k)  < 


i.e.. 


M(k) 


r  i  E(P 
-  E(P;  [-  1,  lj  )  < - 


lLj, 


E(Q(k  +  1;  x);  X(k  +  1))  -  m(k) 

(1  -  X) 

E(P;  [-  1,  1  ]  )  -  m(k) 

(1  -  X ) 


)  -  m(k); 


(1-  X) 

therefore,  lim  M(k)  =  E(P;  [-  1,  l]  ). 

k  — oo  ■* 


The  search  for  the  set  X(k  +  1)  having  obtained  X(k)  is  done  in  the  following  manner:  Start 
at  a  point  x(k,  i)  of  X(k)  and  using  a  step  size  h  we  find  the  direction  of  increasing 
|  Q(k;  x)  -  f(x)  |  and  proceed  in  this  direction.  The  search  continues  until  there  are  two  points 
straddling  the  required  peak  of  the  curve  Q(k;  x)  -  f(x).  The  point  x(k  +  1;  i)  is  then  chosen  as 
the  extreme  of  a  quadratic  interpolation  which  uses  the  last  three  points.  If  the  endpoints -1  and  1 
are  included  inX(k),  we  restrict  the  searchat  these  points  to  one  direction  only,  i.e., to  the  inside 
o£  [-1,  l]. 

There  is  an  important  result  of  a  paper  by  Murnaghan  and  Wrench  "Reference  18)  which 
relates  to  this  search  for  a  set  X(k  +  1).  They  found  that  the  partial  derivative  of  a  coefficient 
of  P(x)  with  respect  to  an  interior  point  of  the  critical  set  corresponding  to  P(x)  is  zero,  and, 
therefore,  this  optimal  critical  set  need  not  be  determined  with  extreme  accuracy.  This  result 
is  valid  only  if  f(x)  is  differentiable  as  well  as  continuous. 
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SECTION  IV 

UPPER  BOUND  FOR  THE  MINIMAX  ERROR 


Let  P(x)  be  the  polynomial  of  degree  less  than  or  equal  to  n  of  best  approximation  to  a  con¬ 
tinuous  function  f(x)  defined  on  j-  1,  1  j . 

Let 

E1  -  max  |  f(x)  -  P(x)  I  for  x  in  |  -  1,  lj  . 

We  will  call  E*  the  minimax  error  for  P(x)  on  £-1,  l| . 

Let  f(x)  have  a  Fourier- Chebycheff  expansion 

oo 

T(x)  =  0.5  aQ  +  Y  ak"^"k 
k=  I 

where  T^(x)  is  the  Chebycheff  polynomial  of  degree  k. 

Theorem  (2):  Let  |a.  j  be  the  subsequence  of  all  non-zero  coefficients  of  T(x).  If  for  all 

k 

k  >n  we  have 


I  a  /a.  I  <  X  <  1 
k  +  1  k 


then  for  k>n  we  have  E '  <  I  a.  1  *  .  ,  where  n  is  the  degree  of  P(x). 

ln  +  1  1  "  X 


oo 


Proof:  If  for  all  k>n  we  have  |  a.  /a.  |  <  X  <  1  then  Y  |  a,  j<oo.  Hence  for  k>n 

\  +1  \  ~  k=l 

00  00 

E*<  I  a. .  T.  (x)  <  Y  I  a,  I, 


]=' 


J  ] 


n+l 


i=' 


) 


n+l 


since  max  T^(x)  =  1  on  1,  1  J. 


But 


00  ?  I 

Y  la  I <  la.  I  (  I  +  X  +  X  +  •  •  )  =  la.  I  — r 

i-'n+l  n+l  n+l 


This  completes  the  proof. 
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SECTION  V 

NUMERICAL  RESULTS 


AH  calculations  were  performed  on  an  IBM  7094  Digital  Computer.  Twenty- seven  bit  accu¬ 
racy  was  used  throughout.  The  linear  programming  code  is  that  of  R.  J.  Clasen  of  the  RAND 
Corporation.  The  subroutine  name  is  RSMSUB,  with  SHARE  Distribution  Number  1281. 

For  most  functions  it  was  found  that  the  sample  independent  variable  set  need  not  contain 
more  than  n+2  points  (where  n  is  the  degree  of  the  approximating  polynomial),  but  in  some 
cases  it  was  found  advantageous  to  use  a  larger  sample  as  this  shortened  the  search  for  the 
critical  set  considerably.  For  all  the  cases  considered  the  end  points  of  the  interval  over  which 
the  approximation  was  obtained  were  included  in  the  critical  set. 


We  can  in  fact  show  that  if  f(x)  has  at  least  (n  +  1)  derivatives  and  if  f^n+^(x)  does  not  change 
signs  in  ^a,  bj ,  then  the  end  points  a  and  b  will  always  be  included  in  the  critical  set. 

Xn+2 

f(x)  is  the  function  we  are  to  approximate,  and  P(x)  is  the  polynomial  of  degree  less  than  or 
equal  to  n  which  best  approximates  f(x)  in  the  interval  [ a,  b).  Let  P(x)  =  Cq+CjX+.  .  .+cnxn. 


be  the  critical  set  of  sample  points.  Let  D(x)  =  f(x)  -P(x),  where 


Let  X  = 


v  x2* 


We  know  that  at  the  interior  points  of  X  the  derivative  of  D(x)  must  be  zero,  i.e., 
D1  (x.)  =  0  i  =  2,  3,  .  .  .  ,  n  +  1 


In  order  to  be  in  the  set  X,  the  points  x^  and  xn  +  2  mus^  either  be  end  points  of  the  interval, 

so  that  x.  =  a  and  x  0  =  b,  or  must  be  such  that  D'(x.)  =  D*(x  ,0  =  0.  We  show  that  given 
1  n+2  '  r  n  +  2'  b 

the  above  conditions  the  latter  cannot  be  true. 


We  know  that  D'(x)  is  equal  to  zero  at  least  n  times  in|  a,  bj.  Repeated  applications  of  Rolle’s 

theorem  show  that  D^n\x)  is  equal  to  zero  at  least  once  in  ^a,  bj.  Since  P(x)  is  of  degree  at 

most  n  we  have  D^n+^(x)  =  f^n+^(x).  Hence,  if  f^n+^(x)  does  not  change  sign  in  [a,  bj,  neither 

will  D'  '(x).  But  then  D^  ^(x)  is  either  monotonic  increasing  or  decreasing,  and  hence  must 

equal  zero  exactly  once.  Then  D^n  ^(x)  will  be  zero  at  most  twice  in  a,  b],  and  continuing  in 

this  way,  we  see  that  D'(x)  will  equal  zero  at  most  n  times,  and,  by  what  was  said  above, 

exactly  n  times.  Then  D'(x)  cannot  be  equal  to  zero  for  x^  and  x^  +  hence,  we  must  have 

x.  =  a  and  x  0  =  b. 

1  n+2 


This  result  is  useful  since  if  f(x)  satisfies  the  required  conditions  on  a,  b  ,  then  we  can 
include  a  and  b  in  our  critical  set  from  the  start. 


Example  1:  We  compute  the  polynomial  P(x)  of  degree  less  than  or  equal  to  3  which  is  the  best 
approximation  to  a  function  f(x)  given  as  a  set  of  ordered  pairs  (x.,  f(x.)).  We  use  as  a  sample 

set  X,  11  equally  spaced  points  ranging  fromO  to  1.  The  values  f(x.)  are  taken  to  be  the  values 

of  the  function  sin(x)  on  the  sample  set  X.  The  Fortran  Library  subroutine  was  used  to  obtain 
the  values  f.  =  f(x.)  =  sin(x.). 
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P(x) 

=  c0  +  CjX  +  c2x2  + 

3 

C3X 

C0 

=  -0.0001472 

c2  =  -0.0193834 

C1 

=  1.00444 

c3  =  -0.143585 

D 

=  max  I  f(x.)  -  P(x.) 
l<i<ll'  1  1 

=  0.0001472 

X. 

1 

f(x.) 

P(*j) 

E(x.)  =  f(x.)  -P(Xj) 

0.0 

0.0 

-0.0001472 

0.0001472 

0.1 

0.0998334 

0.0999592 

-0.0001258 

0.2 

0.198669 

0.198817 

-0.0001472 

0.3 

0.295520 

0.295563 

-0.0000429 

0.4 

0.389418 

0.389338 

0.0000808 

0.5 

0.479426 

0.479278 

0.0001472 

0.6 

0.564642 

0.564524 

0.0001187 

0.7 

0.644218 

0.644213 

0.0000051 

0.8 

0.717356 

0.717483 

-0.0001271 

0.9 

0.783327 

0.784474 

-0.0001472 

1.0 

0.841471 

0.841324 

0.0001472 

Example  2 


Function 


Range 

Approximation  P(x) 
Number  of  Iterations 
E(Q(2,  x);  X(2)) 


x 

\/-/2i t  j  exp  (  -  0.5  t2  )  dt 
-oo 

[4,  4 

c0  +  clx  +  c2x2  +  C3x3 
3 

0.0645903 
0.5 

0.288445 

0.0 

-0.011225 
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The  algorithm  was  terminated  when 

I  E(Q(k  +  1;  x);  X(k  +  1))  -E(Q(k;  x);  X(k))  |  <  0.0000009,  k  =  0,  1 . 

This  took  place  at  k  =  1  for  this  example.  The  set  X(0)  was  as  follows: 
X(0)  =  jx.  :  Xj  =  -4,  x.  +  j  =  x.  +  .4,  i  =  1,  .  .  .  ,  2oj. 

The  critical  set  for  P(x)  was  (±4,  ±2.895,  ±0.92736). 


Example  3 

Function 

Range 

Approximation  P(x) 


Number  of  Iterations 
E(Q(2;  x);  X(2)) 


ln(x  +  1) 

[o,l] 


+ 


Cj  X 


2 


+  C2X 


+  C3X 


+ 


c.x 

4 


+ 


C5X 


3 

0.00000867 

0.00000867 

0.9993 

-0.490747 

0.286714 

-0.133226 

0.0311055 


The  initial  sample  points  were  (0.0,  0.08,  0.28,  0.5,  0.7,  0.9,  1.0).  The  critical  set  for  P(x) 
was  (0.0,  0.060052,  0.234368,  0.484917,  0.73004,  0.926248,  1.0).  The  algorithm  was  terminated 
when  the  difference  between  two  successive  errors  was  less  than  0.0000001. 
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Example  4 
Function 


Range 

Approximation  P(x) 


Number  of  Iterations 
E (Q (2 ;  x);  X(2)) 


f(x)  =  0 

=  0.5 
=  1 


-1  <  x  <  0 
x  =  0 
0  <  x  <  1 


[-i.  i] 

2  3  4  5 

c0  +  cix  +  C2X  +  C3X  +  C4X  +  C5X 


3 

0.5 

0.5 

4.25717 

0.0 

-12.641 

0.0 

9.38384 


The  algorithm  was  started  with  21  equally  spaced  points  in  [-1,  l],  with  the  end  points  in¬ 
cluded.  The  search  for  the  critical  set  required  a  step  size  of  0.0001.  The  critical  set  was 
(±1.0,  ±0.8216,  ±0.3673,  ±0.00000001).  In  this  case  the  maximum  error  took  place  at  eight  points 
rather  than  the  expected  seven.  The  proximity  of  two  of  the  critical  points  makes  the  search 
for  the  critical  set  very  difficult  and  extreme  care  is  required. 
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e  (x  ) 


Figure  3.  Error  Function  for  Example  3 


e  ( x ) 


x 


Figure  4.  Error  Function  for  Example  4 
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APPENDIX 

COMPUTER  PROGRAM 


The  program  was  coded  entirely  in  the  FORTRAN  IV  Programming  Language.  It  consists  of 
four  main  subroutines  and  five  minor  ones.  Their  functions  are  defined  below: 

RELAX  -  serves  as  an  executive  routine  to  the  whole  program  and  develops  certain 
parameters  necessary  for  the  rest  of  the  program. 

APPROX  -  computes  the  polynomial  of  best  approximation  to  the  given  function  over  a 
discrete  point  set. 

SEARCH  -  determines  the  critical  points  and  the  direction  in  which  to  search  for  a  new 
sample  set.  The  search  continues  until  three  points  straddle  the  required  peak 
of  the  error  curve.  These  three  points  ard  the  corresponding  values  of  the  in¬ 
dependent  variable  are  used  in  PEAK. 

PEAK  -  computes  new  critical  point  using  the  Lagrangian  formula  for  a  second  degree 
interpolating  polynomial.  The  three  ordered  pairs  computed  in  SEARCH  are 
used  as  data,  and  the  derivative  of  the  resulting  polynomial  is  solved  for  its 
zero  value  to  obtain  the  required  point. 

SIMPLX  -  used  by  APPROX  to  solve  the  linear  program. 

SETUP  -  used  by  APPROX  to  arrange  the  necessary  input  matrix  to  SIMPLX. 

APPFN  -  used  by  APPROX  to  compute  xj ,  i  =  1 . m 

j  =  0 . n 

FF  -  computes  values  f(x.)  of  given  function. 

ERRFN  -  computes  f(x. )  -P(x. ),  where  P  is  the  approximating  polynomial. 


Limitations  on  the  program  are  (1)  a  maximum  of  201  sample  points  and  (2)  a  maximum  de¬ 
gree  of  9  on  the  approximating  polynomials.  These  limitations  are  easily  extended  if  desired. 

The  calling  sequence  to  the  program  is  as  follows: 

CALL  RELAX  (C,  D,  M,  N,  X,  Y,  ERRJ,  MAXKT,  DELTA,  TOL,  DTOL) 


Output  Parameters: 


C  =  array  of  coefficients  of  desired  polynomial 
D  =  value  of  maximum  error 


ERRJ  =  error  code 
=  0,  no  error 

=  -1,  approximation  not  found  after  allowable  number  of  iterations 
=  -2  error  did  not  increase  on  succeeding  iteration 
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Input  Parameters: 

M  =  number  of  points  in  sample  set 
N  =  n  +  1,  where  n  is  degree  of  polynomial 
X  =  array  of  independent  variable 
Y  =  array  of  dependent  variables  f(x^ ) 

MAXKT  =  maximum  number  of  iterations  allowed 

DELTA  =  step  size  to  be  used  in  search  for  succeeding  critical  sets 

TOL  =  minimum  allowable  value  of  delta 

DTOL  =  routine  terminates  when  absolute  difference  between  succeeding  errors  is 
less  than  DTOL 
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